Nonequilibrium Green's function approach to mesoscopic thermal transport 
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We present a formulation of a nonequilibrium Green's function method for thermal current in 
nanojunction atomic systems with nonlinear interactions. This first-principle approach is applied 
to the calculation of the thermal conductance in carbon nanotube junctions. It is shown that non- 
linearity already becomes important at low temperatures. Nonlinear interactions greatly suppress 
phonon transmission at room temperature. The peak of thermal conductance is found to be around 
400K, in good agreement with experiments. High-order phonon scattering processes are important 
for diffusive heat transport. 
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Thermal transport in materials has been studied for 
a long time, beginning with Joseph Fourier's heat con- 
duction law. However, a microscopic theory is possible 
only after the advent of quantum mechanics JlJ . The ear- 
lier treatments are mostly for bulk systems |2j ■ In recent 
years, motivated by the shrinkage of sizes of electronic 
devices, researchers have paid more attention to the heat 
transport phenomena in meso- and nano-scales Q. Un- 
der such circumstances, some of the concepts have to be 
modified. For example, it has been found that Fourier's 
law is no longer valid for many one-dimensional systems 
0. What to replace it is both interesting theoretically 
and relevant experimentally. 

A number of approaches have been used to study heat 
transport, such as classical molecular dynamics (MD) 
and the Boltzmann-Peierls equation [l| . MD can handle 
nonlincarity, but it is not correct at low temperatures. 
The Boltzmann-Peierls method relies on the concept of 
a phonon distribution function in space which is not any 
more meaningful in nanojunctions where translations! in- 
variance is broken. The Landauer formula takes care 
of the low-temperature limit of ballistic heat transport. 
Some attempts have been made to cover both limits, such 
as a phenomenological investigation using the concept of 
phonon mean free path [f|. Recent works in refs. U 
are efforts from more fundamental points of views, start- 
ing from quantum principles. However, these attempts 
rely on specific models and approximations. Clearly, a 
unified approach valid for the whole temperature range 
is still lacking. 

In this paper we give a theory for heat transport 
in nanojunction using nonequilibrium Green's functions. 
Our approach is an exact, first-principle formulism for 
general models with nonlinear interactions, provided that 
a self-energy can be computed. This technique is used 
extensively in electronic trans por t 0|. Our theory goes 
beyond linear elastic regime P.llfj|| by taking nonlinearity 
perturbatively or through mean-field approximations. A 
comparison of several approximations to the self-energy 
on a one-dimensional (ID) chain suggests that mean-field 
is reliable up to room temperature. We then apply the 



method to short carbon nanotube junctions and compare 
with experimental results. 

We consider an insulating solid where only the vibra- 
tional degrees of freedom are important. The system is 
composed of a left lead and a right lead with an arbi- 
trary junction region. Let the displacement from some 
equilibrium position for the j-th degree of freedom in the 
region a be it", a = L, C, R. The quantum Hamiltonian 
is given by 



H=Y, H a + {u L ) r V 

ct=L,C,R 



z.\t„lc..c + (u c ) T V CR u R + H n , (1) 



where T denotes matrix transpose, H a = ^(u a ) T u a + 
i(ti°) T K a u a , u a is a column vector consisting of all the 
displacement variables in region a, and u a is the corre- 
sponding conjugate momentum. K a is the spring con- 
stant matrix and V LC — (V CL ) T is the coupling matrix 
of the left lead to the central region; similarly for V CR . 
For brevity, we have set all the atomic masses to 1, but 
the formulas are equally applicable to variable masses 
with a transformation 1i i ' X 7 */ Tib n . Also, we'll set the 
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Planck constant fi and Boltzmann constant ks to 1. The 
nonlinear part of the interaction is 



u?ufu%. 



(2) 



Quartic interaction can also be handled. 

A great simplification is possible due to the linear na- 
ture of the leads Hr,, Hr, and the interaction V Ca with 
the central region. Only the central region has nonlinear 
interactions. The leads are assumed semi-infinitely long 
which produce dissipation. A traditional approach for 
many-body problems is to work in second quantization 
framework with the phonon creation and annihilation op- 
erators. Yet, for junction systems without translational 
invariance, we find that the notation and the Green's 
functions will be simpler if we stay in the first quantiza- 
tion and in the coordinate representation [llj . 

The present proposal is parallel to the ideas of the 
nonequilibrium Green's function method for electronic 
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transport. Imagine that at time — oo the system is in 
three separate subsystems in respective thermal equilib- 
rium at the inverse temperature [3 a , a = L,C',R. The 
couplings are switched on adiabatically, so that at time 
t = 0, a steady state is established. A key starting point 
is an expression for the heat current. We begin with the 
definition, 



(3) 



where the decrease in energy of the left lead gives the 
heat flow to the central region. The average is taken 
with respect to an unknown density matrix and will be 
clarified later. By the Heisenberg equation of motion, 
we obtain, at t = 0, I L = {(u L ) T V LC u c ) . The ex- 
pectation value can be expressed in terms of a Green's 
function G£ L (t,t') = ~i(u L (t')u c {t) T ) T . Using the fact 
that operators u and ii are related in Fourier space as 
u[u>] = (—iui)u[u>], we get, 



Il = - 



lj duj. 



(4) 



The next step is to eliminate the reference to the lead 
Green's functions in favor of the Green's functions of the 
central region. We use the contour ordered Green's func- 
tion, defined on a Keldysh contour [12] from — oo to +oo 
and back. The contour ordered Green's function can be 
mapped onto four different normal Green's functions by 
G aa '{t,t') = lim e ^ + G(t+iea,t'+iea'), where a = ±(1), 
and G ++ = G* is the time ordered Green's function, 
G~~ = G* is the anti-time ordered Green's function, 
G+- = G<, and G~ + = G> . The retarded Green's func- 
tion is given by G r — G* — G < , and the advanced by G a = 
G < — G*. These relations also hold for the self energy dis- 
cussed below. In terms of the contour ordered Green's 
function, it can be shown from an equation of motion 
method or direct verification by definition, for our model, 
that G cl (t, t') = J dT"G cc (r, r")V CL g L (r", r% where 
the integral is along the contour. The function is the 
contour ordered Green's function for the semi-infinite free 
left lead in equilibrium at (3l, e.g., the retarded Green's 
function in frequency domain is obtained by the solu- 
tion of [(lo + ir/) 2 - K L ]g r L [u] = I, 7] -»• 0+, where / 
is an identity matrix. Using the Langreth theorem and 
transforming to Fourier space, the above relation gives 
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final expression for the energy current is 
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where the self-energy due to the interaction with the lead 
is Y>l = V gjjV . For notational simplicity, we have 
omitted the subscript G on the Green's functions denot- 
ing the central region. 

Next, we need a method to compute the full Green's 
functions. The perturbative/diagrammatic expansion is 



used to derive Dyson equations. We can treat both the 
coupling V Ca and the nonlinear term H n as perturba- 
tions, or consider only the nonlinearity as a perturbation. 
The latter is simpler in terms of organization. The con- 
tour ordered Green's function is expressed in interaction 
picture: 

G jk (r,r') = -i{T T u?(T)u»(T')) 

= ^(T rU }(r)4(r') e - l / H " (r " )dr ")o, (6) 

where the displacements refer to the central region, the 
operators in the top line are in Heisenberg picture; T T is 
the contour order operator. The average (• • -)o is with re- 
spect to the density matrix of the nonequilibrium steady 
state when H n = 0. Its explicit form is not known, but 
the Wick theorem is still valid. The Green's function Go 
of the linear system can be computed from that of the 
free subsystems: 

g o(t, r')=gc(T, T'yi-JdTidT2gc(r, ti)T,(t 1 ,t 2 )Gq(t 2 ,t'), 

(7) 

where S = El + ^r- The full nonlinear Green's function 
has three types of diagrams in a perturbation expansion. 
Diagrams with loops disconnected from the two terminals 
are zero and can be dropped. There is another class of di- 
agrams where the two terminals are not connected. Such 
diagrams are not zero, but are constants in r. They rep- 
resent a thermal expansion effect, and do not contribute 
to the heat current in Eq. (01 . Finally, the connected 
part of the Green's function satisfies a similar contour 
ordered Dyson equation relating G c to Go using nonlin- 
ear self-energy S n . In ordinary Green's functions and in 
frequency domain (to ar gum ent suppressed), the Dyson 
equations have solutions 
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(10) 
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We note that when there is no nonlinearity, £„ = 0, 
G = Go = G c , the heat current formula, Eq. (JSJ, can be 
simplified to the Landauer formula, 
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duujf [uj] (f L - f R ) , 



(12) 



with the energy transmission coefficient given by the 
Caroli formula [H 0, f[u] = Tr {G T T L G a T R ) , 
T a — i(Sa — where f a is the Bose- 

Einstein distribution function at (3 a . In order 
to facilitate comparison with linear results, we de- 
fine an effective energy transmission by T c ff[a->] = 
iTr{(G r - G«)(E< - £<) +zG<(T R - T L )} /(f L - f R ). 
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FIG. 1: Feynman diagrams for the interaction self-energy E n . 
Each long line corresponds to a propagator Go(t,t')\ each 
vertex is associated with the interaction strength Tijk- All 
internal site indices are summed and contour time variables 
integrated. The number in front of a graph is the factor mul- 
tiplying the graph value. 



The function T c g [ui] is real and even in w. Such effective 
transmission is temperature-dependent. 

So far the equations above are all exact. Several 
ways of approximating the nonlinear contribution to self- 
energy are possible. We can simply truncate the dia- 
grams 15] in the perturbative expansion for the self- 
energy. The diagrams for £„ to 0(T^- k ) are shown in 
Fig. ^ These diagrams are still in the contour variable 
t. For practical calculation, they have to be changed to 
real time t in terms of G CTCr and Fourier transformed to 
the frequency domain. For example, the leading order 
contribution (first two diagrams) to the nonlinear self- 
energy is: 



n,jk 
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Mean-field theory can be obtained by replacing Go by G, 
and the equations are solved iteratively. 

A general program is implemented based on Eq. J^J, 
Eq. JBJ to fTTJl . and Eq. fl"3|) . The surface Green's func- 
tion g r L is computed using an efficient recursive method 
[l6|. In numerical calculation of the Green's functions, it 
is important to keep a small but finite r], as the functions 
are singular in the limit 77 — ^ + . In addition, on-site 
potentials are applied to the leads to make the system 
stable. 

We first test various approximations on a ID junction 
with parameters comparable to that of a carbon chain. 
The system consists of harmonic leads and a junction 
part with harmonic interactions plus cubic interactions 
of the form (l/3)i^(uj ~ u j+i) 3 °f Fermi-Pasta-Ulam 
type. Figure [3 presents the results for a 3-atom junction 
system. We discuss the effect of nonlinearity to thermal 
transport. As we can see, adding nonlinear contribu- 
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FIG. 2: Thermal conductance as a function of temperature 
for a ID junction with three atoms. The harmonic spring con- 
stants are k L = 1.56, k R = 1.44, k c = 1.38 (eV/(A 2 amu)). 
The nonlinear strength is t = 1.8 eV/(A 3 amu 3 ^ 2 ). Small on- 
site quadratic potentials are applied to the leads with spring 
constants kT site = 0.01, and it™ site = 0.02 (eV/(A 2 amu)). 



tions suppresses the thermal conductance at high tem- 
peratures. The deviation from ballistic transport starts 
around 200 Kelvin. As the temperature rises further, we 
expect that the higher order graphs become important. 
The high order calculations are rather expensive, with 
computational complexity scaled as 0(N 4 M 2 ) where N 
is system size and M is sampling points in frequency. To 
partially take into account the higher order contributions 
but still keep the computation within reasonable limit for 
large systems, we find a mean-field theory is most satis- 
factory. In this version of mean-field treatment, we con- 
sider the leading diagrams of 0(T 2 ), and replace Go by G 
only for the first diagram. The equations are then solved 
iteratively. Good agreement with 0(T A ) result is found 
for temperatures up to 400 K for the ID chain. Thus we 
expect that the mean-field theory can give excellent re- 
sults for moderately high temperatures. However, both 
the perturbative results and mean-field one break down 
at sufficiently high temperature, as the cubic nonlinear- 
ity is only metastable. To predict correctly diffusive be- 
havior at high temperatures, the quartic interaction and 
higher order graphs are essential. 

We now turn to the carbon nanotube junctions. The 
leads and the junction are of the same diameter (8,0) 
nanotube. The only difference is that the junction part 
has cubic nonlinear interactions while the leads are per- 
fectly harmonic. The values of T^k are derived from 
the Brenner potential by finite differences. For compu- 
tational efficiency, small values of T^-fe are truncated to 
zero. Figure shows the ballistic transmission coefficient 
when the nonlinearity is set to zero, and is compared 
with the effective transmission due to the leading non- 
linear effects (the first two graphs). It can be seen that 
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Ballistic 

Effective Transmission 




can estimate that the thermal conductivity at the exper- 
imental accessible length is about 2000 W/(mK). This is 
qualitatively in agreement with experimental values |l7| . 

We proposed a fully quantum mechanical approach for 
computing heat current of solid junctions with nonlinear 
interactions. We have demonstrated the method with 
ID and molecular junction systems. The perturbation 
expansion for self-energy works well up to room temper- 
atures. However, it is still a challenge to find efficient 
and good approximation for the self-energy at high tem- 
peratures. 

This work is supported in part by a Faculty Research 
Grant of the National University of Singapore. 
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FIG. 3: The ballistic transmission and the effective trans- 
mission at 300K for an (8,0) one-unit-cell carbon nanotube 
junction. 
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FIG. 4: Thermal conductance of (8,0) carbon nanotube junc- 
tion with one unit cell (0.426 nm). The inset shows the ther- 
mal conductivity as a function of tube length. 

the nonlinearity greatly suppresses the thermal transmis- 
sion, particularly at low frequencies. The discontinuity is 
smeared out to more smooth curve due to thermal effect. 

The temperature and length dependence of the nan- 
otube thermal conductance (conductivity) is shown in 
Fig. 2| The mean-field theory result gives a peak in the 
conductance around 400K. This behavior agrees with ex- 
periments 17) and is in contrast with MD results which 
tend to give peaks at much lower temperatures 01 ■ The 
inset shows the thermal conductivity calculated from 
the conductance. The cross-section is defined as 7r<i 2 /4, 
where d is the diameter of the tube. If we assume that 
the transport up to 1 /jm is still quasi-ballistic, then we 
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